1 확진자 데이터
COVID 팩키지에서 한국에서 코로나19 검사결과 확진된 시계열 데이터를 얻는다. 시계열 데이터 예측을 위한 기초 데이터를 준비한다.
library(tidyverse)
library(COVID19)
library(timetk)
korea_dat <- covid19("KOR", level = 1)We have invested a lot of time and effort in creating COVID-19 Data Hub, please cite the following when using it:
Guidotti, E., Ardia, D., (2020), "COVID-19 Data Hub", Journal of Open
Source Software 5(51):2376, doi: 10.21105/joss.02376.
A BibTeX entry for LaTeX users is
@Article{,
title = {COVID-19 Data Hub},
year = {2020},
doi = {10.21105/joss.02376},
author = {Emanuele Guidotti and David Ardia},
journal = {Journal of Open Source Software},
volume = {5},
number = {51},
pages = {2376},
}
To retrieve citation and metadata of the data sources see ?covid19cite. To hide this message use 'verbose = FALSE'.korea_df <- korea_dat %>%
ungroup() %>%
dplyr::select(날짜=date, 누적검사자=tests, 누적확진자=confirmed, 누적회복자=recovered, 누적사망자=deaths) %>%
timetk::pad_by_time(.date_var = 날짜, .by = "day", .pad_value = NA) %>%
mutate(누적검사자 = timetk::ts_impute_vec(누적검사자, period = 1) %>% ceiling(.),
누적회복자 = ifelse(is.na(누적회복자), 0, 누적회복자),
누적사망자 = ifelse(is.na(누적사망자), 0, 누적사망자))
korea_daily_df <- korea_df %>%
mutate(검사자 = 누적검사자 - lag(누적검사자, n = 1L),
확진자 = 누적확진자 - lag(누적확진자, n = 1L),
회복자 = 누적회복자 - lag(누적회복자, n = 1L),
사망자 = 누적사망자 - lag(누적사망자, n = 1L) ) %>%
drop_na() %>%
select(-contains("누적"))
confirmed_daily_tbl <- korea_daily_df %>%
select(날짜, 확진자) %>%
timetk::pad_by_time(.date_var = 날짜, .by="day", .pad_value = NA)
confirmed_daily_tbl # A tibble: 339 x 2
날짜 확진자
<date> <int>
1 2020-01-23 0
2 2020-01-24 1
3 2020-01-25 0
4 2020-01-26 1
5 2020-01-27 1
6 2020-01-28 0
7 2020-01-29 0
8 2020-01-30 0
9 2020-01-31 7
10 2020-02-01 1
# ... with 329 more rows2 시각화
timetk 팩키지 plot_time_series() 함수를 사용해서 확진자 추세를 시각화한다.
confirmed_daily_tbl %>%
plot_time_series(.date_var = 날짜,
.value = 확진자,
.line_color = "#ff0000",
.title = "코로나19 일별 확진자 추세")3 모형 데이터
3.1 전체 데이터
confirmed_daily_tbl 데이터를 history_tbl 데이터프레임으로 저장시켜 명확히 의미를 부여한다. timetk 팩키지 future_frame() 함수로 예측할 기간(“1달”)으로 데이터프레임을 확장하여 관측된 시계열과 예측 시계열을 결합시킨 전체 데이터프레임 데이터를 준비한다. forecast_tbl 로 예측할 데이터프레임을 따로 떼어둔다.
# 전체 시계열 데이터 ----
full_tbl <- confirmed_daily_tbl %>%
timetk::future_frame(.length_out = "1 month",
.bind_data = TRUE)
# 관측 시계열 데이터 ----
history_tbl <- full_tbl %>%
filter(!is.na(확진자))
# 예측 시계열 데이터 ----
forecast_tbl <- full_tbl %>%
filter(is.na(확진자))
full_tbl %>%
plot_time_series(.date_var = 날짜,
.value = 확진자,
.title = "코로나19 일별 확진자수")3.2 훈련/시험 데이터
time_series_split() 함수를 사용해서 훈련/시험 표본으로 데이터를 구분한다. assess = 28을 지정하여 향후 28일 후를 예측하는 모형을 개발한다.
splits <- history_tbl %>%
time_series_split(날짜, assess = 28, cumulative = TRUE)
# time_series_split(날짜, initial = "3 month", assess = 28, cumulative = FALSE)
splits<Analysis/Assess/Total>
<311/28/339>상기 훈련/시험 데이터 구분한 후 시각화를 통해 확인해보자. timetk::plot_time_series_cv_plan() 함수를 사용한다.
splits %>%
tk_time_series_cv_plan() %>%
timetk::plot_time_series_cv_plan(날짜, 확진자)4 피처 공학
tidymodels 생태계의 recipe 팩키지 recipe() 함수를 사용해서 피쳐 공학(feature engineering) 작업을 수행한다. 추후, 다양한 피처를 작업할 예정이라 대략적인 틀만 잡아둔다.
library(tidymodels)
recipe_spec <- recipes::recipe(확진자 ~ ., data = training(splits))
recipe_spec %>% prep() %>% juice()# A tibble: 311 x 2
날짜 확진자
<date> <int>
1 2020-01-23 0
2 2020-01-24 1
3 2020-01-25 0
4 2020-01-26 1
5 2020-01-27 1
6 2020-01-28 0
7 2020-01-29 0
8 2020-01-30 0
9 2020-01-31 7
10 2020-02-01 1
# ... with 301 more rows5 모형적합
parsnip 팩키지의 linear_reg() 를 통해 선형회귀(?) 모형을 시계열 데이터에 날짜를 독립변수로 넣어 말도 되지 않지만 가장 단순하게 예측 모형을 적합시킨다. 이를 위해서 workflow() 를 사용하는데 기본적으로 피처 공학의 recipe_spec이 필요하고 model_spec의 모형이 필요하다.
그리고 fit() 함수를 사용해서 적합을 시킨다.
model_spec <- linear_reg(
mode = "regression") %>%
set_engine("lm")
wkfl_fit_lm <- workflow() %>%
add_recipe(recipe_spec) %>%
add_model(model_spec) %>%
fit(training(splits))
wkfl_fit_lm== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: linear_reg()
-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps
-- Model -----------------------------------------------------------------------
Call:
stats::lm(formula = ..y ~ ., data = data)
Coefficients:
(Intercept) 날짜
-5725.9014 0.3164 6 모형 평가
workflow 객체를 modeltime 팩키지 modeltime_table() 함수에 넣어 객체로 만든 후에 modeltime_accuracy() 함수를 사용해서 모형 성능을 파악한다.
library(modeltime)
model_tbl <- modeltime_table(
wkfl_fit_lm
)
model_tbl %>%
modeltime::modeltime_accuracy(testing(splits))# A tibble: 1 x 9
.model_id .model_desc .type mae mape mase smape rmse rsq
<int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 LM Test 661. 78.4 6.48 130. 701. 0.8026.1 시각화
시험데이터를 통해 모형으로 예측한 값을 시각화한다. \(R^2\) 값도 80%가 나오는데 시각화를 하니 시험데이터에 대한 적합이 데이터와 동떨어진 것을 확인할 수 있다.
calibration_tbl <- model_tbl %>%
modeltime_calibrate(
new_data = testing(splits)
)
calibration_tbl %>%
modeltime_forecast(
new_data = testing(splits),
actual_data = history_tbl,
conf_interval = 0.10
) %>%
plot_modeltime_forecast(
.legend_max_width = 60,
.legend_show = FALSE,
.conf_interval_show = TRUE,
.conf_interval_alpha = 0.20,
.conf_interval_fill = "lightblue",
.title = "코로나19 확진자 1개월 예측 - 선형회귀"
)7 확진자 예측
앞선 모형은 history_tbl 을 훈련/시험 데이터에 대해 적합을 시킨 것이라 … 이제 시간을 확대하여 modeltime_refit() 함수를 사용해서 모형을 전체 데이터에 대해 다시 적합시킨다.
refit_tbl <- calibration_tbl %>%
modeltime_refit(data = history_tbl)마지막으로 앞서 구축된 모형을 바탕으로 현재까지 입수된 데이터를 바탕으로 한달 후를 예측한다.
refit_tbl %>%
modeltime_forecast(
new_data = forecast_tbl,
actual_data = history_tbl,
conf_interval = 0.5
) %>%
plot_modeltime_forecast(
.legend_max_width = 25,
.conf_interval_fill = "lightblue",
.conf_interval_alpha = 0.7,
.interactive = FALSE
)
8 모형과 데이터 저장
마지막 단계로 데이터와 모형을 저장하여 다음 단계로 나가기 위한 작업을 수행한다.
# 데이터 ----
full_tbl %>%
write_rds("data/full_tbl.rds")
# 모형 ----
wkfl_fit_lm %>%
write_rds("data/wkfl_fit_lm.rds")데이터 과학자 이광춘 저작
kwangchun.lee.7@gmail.com